clear all;
clc;
format long;
%% parameter settings 
rng('shuffle');

NP = 800;        % population size 
srmin = 0.05;
srmax = 0.5;
run_times = 30;  % independent experiment times
Dim = 30;        % dimension
lb = -100;       % lower bound
ub = 100;        % upper bound
pop_exp = 2;     % cs decrease rate
top_exp = 0.1;   % sr decrease rate

FEs = Dim*1e4;
lb_matrix = lb(ones(1,NP),ones(1,Dim));
ub_matrix = ub(ones(1,NP),ones(1,Dim));
top_ub = ub(ones(1,Dim));
top_lb = lb(ones(1,Dim));
local_search = 2; 
sigma = repmat(1e-4,1,Dim);
folder = './results';
mkdir(folder);

for func_id = 1:30    % 30 benchmark functions
%% store data
optima = func_id*100; % optimum
file_str=[num2str(func_id),'_fitness_EDA','.txt'];
path=['./results/',file_str];
fileID=fopen(path,'w');

    for run = 1:run_times
    %% initail population
        fit_counter = 0;
        record = 2;
        best_list = [];
        ini_pop = lb + (ub-lb)*rand(NP,Dim);
        fit = cec14_func(ini_pop',func_id)';
        fit_counter = fit_counter + NP;
        [~,sorts] = sort(fit);
        pop = ini_pop(sorts(1:NP),:);
        fit = fit(sorts(1:NP));
        
        best_fit = fit(1);
        best_ind = pop(1,:);
        best_list = best_fit;
        par_pop = pop;
        par_fit = fit;
    %% the main loop
        while(fit_counter < FEs)
            if(fit_counter > (record - 1) * Dim*100)
                best_list = [best_list best_fit];
                record = record + 1;
                disp(best_fit);
            end
            
            s=ceil(NP*srmax)-ceil((NP*srmax-NP*srmin)*power((fit_counter/FEs),top_exp));
            sc = NP-ceil((NP-NP*srmin)*power((fit_counter/FEs),pop_exp)); 
            u = mean(pop(1:s,:),1);  % calculate mean vector
            C = ((pop(1: sc,:)-u(ones(1,sc),:))'*(pop(1: sc,:)-u(ones(1,sc),:))) / (sc-1); % calculate covariance matrix

            C = triu(C) + transpose(triu(C,1));
            [B,D] = eig(C);
            if(max(diag(D)) > 1e15*min(diag(D)))          
                tmp = max(diag(D)) / 1e15 - min(diag(D)); 
                C = C + tmp * eye(Dim);
                [B,D] = eig(C);
            end
            D = diag(sqrt(diag(D)));
            pop = u(ones(1,NP),:) + randn(NP,Dim)*(D'* B');     % sample new population

            index = pop > ub_matrix;                                 % bound control 
            pop(index) = 2*ub_matrix(index) - pop(index);
            index = pop < lb_matrix;
            pop(index) = 2*lb_matrix(index) - pop(index);

            fit = cec14_func(pop',func_id)';       
            fit_counter = fit_counter + NP;

            tmp_par_pop = pop;  %store father pop
            tmp_par_fit = fit;

            tmp_pop = par_pop;
            tmp_pop(NP+1:2*NP,:) = pop;
            tmp_fit = par_fit;
            tmp_fit(NP+1:2*NP) = fit;

            par_pop = tmp_par_pop;
            par_fit = tmp_par_fit;

            [~,sorts] = sort(tmp_fit);
            pop = tmp_pop(sorts(1:NP),:); % 2np--np 
            fit = tmp_fit(sorts(1:NP));

            if(best_fit < fit(1))         % store global best fitness
                pop(1,:) = best_ind;
            else
                best_ind = pop(1,:);
                best_fit = fit(1);
            end 
       %% local search
            for i = 1:local_search
                local_ind = normrnd(best_ind,sigma);
                
                index = local_ind > top_ub;                               % bound control 
                local_ind(index) = 2*top_ub(index) - local_ind(index);
                index = local_ind <  top_lb;
                local_ind(index) = 2*top_lb(index) - local_ind(index);
                
                local_fit = cec14_func(local_ind',func_id);
                if local_fit < best_fit
                    best_fit = local_fit;
                    best_ind = local_ind;
                    pop(1,:) = local_ind;
                end
            end
            fit_counter = fit_counter + local_search; 
        end
        best_list = [best_list best_fit];
        fprintf(fileID,'%e\t',best_list-optima);    % save data
        fprintf(fileID,'\r\n'); 
    end
    disp(func_id);
    fclose(fileID);
end